

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>plfit.plfit_py &mdash; agpy 0.1.2 documentation</title>
    
    <link rel="stylesheet" href="../../_static/extra.css" type="text/css" />
    <link rel="stylesheet" href="../../_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    '../../',
        VERSION:     '0.1.2',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="../../_static/jquery.js"></script>
    <script type="text/javascript" src="../../_static/underscore.js"></script>
    <script type="text/javascript" src="../../_static/doctools.js"></script>
    <link rel="top" title="agpy 0.1.2 documentation" href="../../index.html" />
    <link rel="up" title="Module code" href="../index.html" />
     
    <script type="text/javascript">

      var _gaq = _gaq || [];
      _gaq.push(['_setDomainName', 'pyspeckit.bitbucket.org']);
      _gaq.push(['_setAllowHash', false]);
      _gaq.push(['_trackPageview']);


    </script>
    <link rel="stylesheet" type="text/css" href="../../_static/extra.css" />
  </head>
  <body>
    <div class="header-wrapper">
      <div class="header">
        <h1><a href="../../index.html">agpy 0.1.2 documentation</a></h1>
        <div class="rel">
          <a href="http://agpy.googlecode.com">agpy Home </a>  |
          <a href=../../index.html>Docs Home </a>  |
          <a href="http://code.google.com/p/agpy/w/list">Wiki</a>  |
          <a href=../../search.html>Search </a>
        </div>
       </div>
    </div>

    <div class="content-wrapper">
      <div class="content">
        <div class="document">
            
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body">
            
  <h1>Source code for plfit.plfit_py</h1><div class="highlight"><pre>
<span class="c"># intended to implement a power-law fitting routine as specified in.....</span>
<span class="c"># http://www.santafe.edu/~aaronc/powerlaws/</span>
<span class="c">#</span>
<span class="c"># The MLE for the power-law alpha is very easy to derive given knowledge</span>
<span class="c"># of the lowest value at which a power law holds, but that point is </span>
<span class="c"># difficult to derive and must be acquired iteratively.</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">plfit.py - a python power-law fitter based on code by Aaron Clauset</span>
<span class="sd">http://www.santafe.edu/~aaronc/powerlaws/</span>
<span class="sd">http://arxiv.org/abs/0706.1062 &quot;Power-law distributions in empirical data&quot; </span>
<span class="sd">Requires pylab (matplotlib), which requires numpy</span>

<span class="sd">example use:</span>
<span class="sd">from plfit import plfit</span>

<span class="sd">MyPL = plfit(mydata)</span>
<span class="sd">MyPL.plotpdf(log=True)</span>


<span class="sd">&quot;&quot;&quot;</span>

<span class="kn">import</span> <span class="nn">time</span>
<span class="kn">import</span> <span class="nn">random</span>
<span class="kn">import</span> <span class="nn">math</span>

<div class="viewcode-block" id="plfit"><a class="viewcode-back" href="../../plfit.html#plfit.plfit_py.plfit">[docs]</a><span class="k">class</span> <span class="nc">plfit</span><span class="p">:</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    A Python implementation of the Matlab code http://www.santafe.edu/~aaronc/powerlaws/plfit.m</span>
<span class="sd">    from http://www.santafe.edu/~aaronc/powerlaws/</span>

<span class="sd">    See A. Clauset, C.R. Shalizi, and M.E.J. Newman, &quot;Power-law distributions</span>
<span class="sd">    in empirical data&quot; SIAM Review, 51, 661-703 (2009). (arXiv:0706.1062)</span>
<span class="sd">    http://arxiv.org/abs/0706.1062</span>

<span class="sd">    The output &quot;alpha&quot; is defined such that p(x) ~ (x/xmin)^-alpha</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span><span class="n">x</span><span class="p">,</span><span class="o">**</span><span class="n">kwargs</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Initializes and fits the power law.  Can pass &quot;quiet&quot; to turn off </span>
<span class="sd">        output (except for warnings; &quot;silent&quot; turns off warnings)</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="n">neg</span> <span class="o">=</span> <span class="p">[</span><span class="n">i</span><span class="o">&lt;</span><span class="mi">0</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">x</span><span class="p">]</span>
        <span class="k">if</span> <span class="nb">any</span><span class="p">(</span><span class="n">neg</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="k">print</span> <span class="s">&quot;Removed </span><span class="si">%i</span><span class="s"> negative points&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="nb">sum</span><span class="p">(</span><span class="n">neg</span><span class="p">))</span>
            <span class="n">x</span> <span class="o">=</span> <span class="p">[</span><span class="n">i</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">x</span> <span class="k">if</span> <span class="n">i</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">]</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">x</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">plfit</span><span class="p">(</span><span class="o">**</span><span class="n">kwargs</span><span class="p">)</span>


<div class="viewcode-block" id="plfit.alpha_"><a class="viewcode-back" href="../../plfit.html#plfit.plfit_py.plfit.alpha_">[docs]</a>    <span class="k">def</span> <span class="nf">alpha_</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span><span class="n">x</span><span class="p">):</span>
        <span class="k">def</span> <span class="nf">alpha</span><span class="p">(</span><span class="n">xmin</span><span class="p">,</span><span class="n">x</span><span class="o">=</span><span class="n">x</span><span class="p">):</span>
            <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">            given a sorted data set and a minimum, returns power law MLE fit</span>
<span class="sd">            data is passed as a keyword parameter so that it can be vectorized</span>
<span class="sd">            &quot;&quot;&quot;</span>
            <span class="n">x</span> <span class="o">=</span> <span class="p">[</span><span class="n">i</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">x</span> <span class="k">if</span> <span class="n">i</span><span class="o">&gt;=</span><span class="n">xmin</span><span class="p">]</span>
            <span class="n">n</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
            <span class="n">divsum</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">([</span><span class="n">math</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">i</span><span class="o">/</span><span class="n">xmin</span><span class="p">)</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">x</span><span class="p">])</span>
            <span class="k">if</span> <span class="n">divsum</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span> <span class="k">return</span> <span class="nb">float</span><span class="p">(</span><span class="s">&#39;inf&#39;</span><span class="p">)</span>
            <span class="n">a</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">n</span><span class="p">)</span> <span class="o">/</span> <span class="n">divsum</span>
            <span class="k">return</span> <span class="n">a</span>
        <span class="k">return</span> <span class="n">alpha</span>
</div>
<div class="viewcode-block" id="plfit.kstest_"><a class="viewcode-back" href="../../plfit.html#plfit.plfit_py.plfit.kstest_">[docs]</a>    <span class="k">def</span> <span class="nf">kstest_</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span><span class="n">x</span><span class="p">):</span>
        <span class="k">def</span> <span class="nf">kstest</span><span class="p">(</span><span class="n">xmin</span><span class="p">,</span><span class="n">x</span><span class="o">=</span><span class="n">x</span><span class="p">):</span>
            <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">            given a sorted data set and a minimum, returns power law MLE ks-test w/data</span>
<span class="sd">            data is passed as a keyword parameter so that it can be vectorized</span>

<span class="sd">            The returned value is the &quot;D&quot; parameter in the ks test...</span>
<span class="sd">            &quot;&quot;&quot;</span>
            <span class="n">x</span> <span class="o">=</span> <span class="p">[</span><span class="n">i</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">x</span> <span class="k">if</span> <span class="n">i</span><span class="o">&gt;=</span><span class="n">xmin</span><span class="p">]</span>
            <span class="n">n</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">))</span>
            <span class="k">if</span> <span class="n">n</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span> <span class="k">return</span> <span class="nb">float</span><span class="p">(</span><span class="s">&#39;inf&#39;</span><span class="p">)</span>
            <span class="n">divsum</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">([</span><span class="n">math</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">i</span><span class="o">/</span><span class="n">xmin</span><span class="p">)</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">x</span><span class="p">])</span>
            <span class="k">if</span> <span class="n">divsum</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span> <span class="k">return</span> <span class="nb">float</span><span class="p">(</span><span class="s">&#39;inf&#39;</span><span class="p">)</span>
            <span class="n">a</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">n</span><span class="p">)</span> <span class="o">/</span> <span class="n">divsum</span>
            <span class="n">cx</span> <span class="o">=</span> <span class="p">[</span><span class="nb">float</span><span class="p">(</span><span class="n">i</span><span class="p">)</span><span class="o">/</span><span class="nb">float</span><span class="p">(</span><span class="n">n</span><span class="p">)</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">n</span><span class="p">)]</span>
            <span class="n">cf</span> <span class="o">=</span> <span class="p">[</span><span class="mi">1</span><span class="o">-</span><span class="p">(</span><span class="n">xmin</span><span class="o">/</span><span class="n">i</span><span class="p">)</span><span class="o">**</span><span class="n">a</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">x</span><span class="p">]</span>
            <span class="n">ks</span> <span class="o">=</span> <span class="nb">max</span><span class="p">([</span><span class="nb">abs</span><span class="p">(</span><span class="n">a</span><span class="o">-</span><span class="n">b</span><span class="p">)</span> <span class="k">for</span> <span class="n">a</span><span class="p">,</span><span class="n">b</span> <span class="ow">in</span> <span class="nb">zip</span><span class="p">(</span><span class="n">cf</span><span class="p">,</span><span class="n">cx</span><span class="p">)])</span>
            <span class="k">return</span> <span class="n">ks</span>
        <span class="k">return</span> <span class="n">kstest</span>
    
</div>
<div class="viewcode-block" id="plfit.plfit"><a class="viewcode-back" href="../../plfit.html#plfit.plfit_py.plfit.plfit">[docs]</a>    <span class="k">def</span> <span class="nf">plfit</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span><span class="n">nosmall</span><span class="o">=</span><span class="bp">True</span><span class="p">,</span><span class="n">finite</span><span class="o">=</span><span class="bp">False</span><span class="p">,</span><span class="n">quiet</span><span class="o">=</span><span class="bp">False</span><span class="p">,</span><span class="n">silent</span><span class="o">=</span><span class="bp">False</span><span class="p">,</span><span class="n">usefortran</span><span class="o">=</span><span class="bp">False</span><span class="p">,</span><span class="n">usecy</span><span class="o">=</span><span class="bp">False</span><span class="p">,</span>
            <span class="n">xmin</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">verbose</span><span class="o">=</span><span class="bp">False</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        A Python implementation of the Matlab code http://www.santafe.edu/~aaronc/powerlaws/plfit.m</span>
<span class="sd">        from http://www.santafe.edu/~aaronc/powerlaws/</span>

<span class="sd">        See A. Clauset, C.R. Shalizi, and M.E.J. Newman, &quot;Power-law distributions</span>
<span class="sd">        in empirical data&quot; SIAM Review, 51, 661-703 (2009). (arXiv:0706.1062)</span>
<span class="sd">        http://arxiv.org/abs/0706.1062</span>

<span class="sd">        nosmall is on by default; it rejects low s/n points</span>
<span class="sd">        can specify xmin to skip xmin estimation</span>

<span class="sd">        There are 3 implementations of xmin estimation.  The fortran version is fastest, the C (cython)</span>
<span class="sd">        version is ~10% slower, and the python version is ~3x slower than the fortran version.</span>
<span class="sd">        Also, the cython code suffers ~2% numerical error relative to the fortran and python for unknown</span>
<span class="sd">        reasons.</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="n">x</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">data</span>
        <span class="n">z</span> <span class="o">=</span> <span class="nb">sorted</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
        <span class="n">t</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span>
        <span class="n">xmins</span> <span class="o">=</span> <span class="nb">sorted</span><span class="p">(</span><span class="nb">set</span><span class="p">(</span><span class="n">z</span><span class="p">))</span>
        <span class="n">argxmins</span> <span class="o">=</span> <span class="p">[</span><span class="n">z</span><span class="o">.</span><span class="n">index</span><span class="p">(</span><span class="n">i</span><span class="p">)</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">xmins</span><span class="p">]</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">_nunique</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">xmins</span><span class="p">)</span>
        <span class="k">if</span> <span class="n">xmin</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
            <span class="n">av</span>  <span class="o">=</span> <span class="nb">map</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">alpha_</span><span class="p">(</span><span class="n">z</span><span class="p">),</span><span class="n">xmins</span><span class="p">)</span>
            <span class="n">dat</span> <span class="o">=</span> <span class="nb">map</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">kstest_</span><span class="p">(</span><span class="n">z</span><span class="p">),</span><span class="n">xmins</span><span class="p">)</span>
            <span class="n">sigma</span> <span class="o">=</span> <span class="p">[(</span><span class="n">a</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span><span class="o">/</span><span class="n">math</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">z</span><span class="p">)</span><span class="o">-</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span> <span class="k">for</span> <span class="n">a</span><span class="p">,</span><span class="n">i</span> <span class="ow">in</span> <span class="nb">zip</span><span class="p">(</span><span class="n">av</span><span class="p">,</span><span class="n">argxmins</span><span class="p">)]</span>
            <span class="k">if</span> <span class="n">nosmall</span><span class="p">:</span>
                <span class="c"># test to make sure the number of data points is high enough</span>
                <span class="c"># to provide a reasonable s/n on the computed alpha</span>
                <span class="n">goodvals</span> <span class="o">=</span> <span class="p">[</span><span class="n">s</span><span class="o">&lt;</span><span class="mf">0.1</span> <span class="k">for</span> <span class="n">s</span> <span class="ow">in</span> <span class="n">sigma</span><span class="p">]</span>
                <span class="k">if</span> <span class="bp">False</span> <span class="ow">in</span> <span class="n">goodvals</span><span class="p">:</span> 
                    <span class="n">nmax</span> <span class="o">=</span> <span class="n">goodvals</span><span class="o">.</span><span class="n">index</span><span class="p">(</span><span class="bp">False</span><span class="p">)</span>
                    <span class="n">dat</span> <span class="o">=</span> <span class="n">dat</span><span class="p">[:</span><span class="n">nmax</span><span class="p">]</span>
                    <span class="n">xmins</span> <span class="o">=</span> <span class="n">xmins</span><span class="p">[:</span><span class="n">nmax</span><span class="p">]</span>
                    <span class="n">av</span> <span class="o">=</span> <span class="n">av</span><span class="p">[:</span><span class="n">nmax</span><span class="p">]</span>
                <span class="k">else</span><span class="p">:</span>
                    <span class="k">print</span> <span class="s">&quot;Not enough data left after flagging - using all positive data.&quot;</span>
            <span class="k">if</span> <span class="ow">not</span> <span class="n">quiet</span><span class="p">:</span> <span class="k">print</span> <span class="s">&quot;PYTHON plfit executed in </span><span class="si">%f</span><span class="s"> seconds&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span><span class="o">-</span><span class="n">t</span><span class="p">)</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">_av</span> <span class="o">=</span> <span class="n">av</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">_xmin_kstest</span> <span class="o">=</span> <span class="n">dat</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">_sigma</span> <span class="o">=</span> <span class="n">sigma</span>
            <span class="n">xmin</span>  <span class="o">=</span> <span class="n">xmins</span><span class="p">[</span><span class="n">dat</span><span class="o">.</span><span class="n">index</span><span class="p">(</span><span class="nb">min</span><span class="p">(</span><span class="n">dat</span><span class="p">))]</span> 
        <span class="n">z</span>     <span class="o">=</span> <span class="p">[</span><span class="n">i</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">z</span> <span class="k">if</span> <span class="n">i</span> <span class="o">&gt;=</span> <span class="n">xmin</span><span class="p">]</span>
        <span class="n">n</span>     <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">z</span><span class="p">)</span>
        <span class="n">alpha</span> <span class="o">=</span> <span class="mi">1</span> <span class="o">+</span> <span class="n">n</span> <span class="o">/</span> <span class="nb">sum</span><span class="p">(</span> <span class="p">[</span><span class="n">math</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">a</span><span class="o">/</span><span class="n">xmin</span><span class="p">)</span> <span class="k">for</span> <span class="n">a</span> <span class="ow">in</span> <span class="n">z</span><span class="p">]</span> <span class="p">)</span> 
        <span class="k">if</span> <span class="n">finite</span><span class="p">:</span>
            <span class="n">alpha</span> <span class="o">=</span> <span class="n">alpha</span><span class="o">*</span><span class="p">(</span><span class="n">n</span><span class="o">-</span><span class="mf">1.</span><span class="p">)</span><span class="o">/</span><span class="n">n</span><span class="o">+</span><span class="mf">1.</span><span class="o">/</span><span class="n">n</span>
        <span class="k">if</span> <span class="n">n</span> <span class="o">==</span> <span class="mi">1</span> <span class="ow">and</span> <span class="ow">not</span> <span class="n">silent</span><span class="p">:</span>
            <span class="k">print</span> <span class="s">&quot;Failure: only 1 point kept.  Probably not a power-law distribution.&quot;</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">_alpha</span> <span class="o">=</span> <span class="mi">0</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">_alphaerr</span> <span class="o">=</span> <span class="mi">0</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">_likelihood</span> <span class="o">=</span> <span class="mi">0</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">_ks</span> <span class="o">=</span> <span class="mi">0</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">_ks_prob</span> <span class="o">=</span> <span class="mi">0</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">_xmin</span> <span class="o">=</span> <span class="n">xmin</span>
            <span class="k">return</span> <span class="n">xmin</span><span class="p">,</span><span class="mi">0</span>
        <span class="k">if</span> <span class="n">n</span> <span class="o">&lt;</span> <span class="mi">50</span> <span class="ow">and</span> <span class="ow">not</span> <span class="n">finite</span> <span class="ow">and</span> <span class="ow">not</span> <span class="n">silent</span><span class="p">:</span>
            <span class="k">print</span> <span class="s">&#39;(PLFIT) Warning: finite-size bias may be present. n=</span><span class="si">%i</span><span class="s">&#39;</span> <span class="o">%</span> <span class="n">n</span>
        <span class="c"># ks = max(abs( numpy.arange(n)/float(n) - (1-(xmin/z)**(alpha-1)) ))</span>
        <span class="n">ks</span> <span class="o">=</span> <span class="nb">max</span><span class="p">(</span> <span class="p">[</span><span class="nb">abs</span><span class="p">(</span> <span class="n">i</span><span class="o">/</span><span class="nb">float</span><span class="p">(</span><span class="n">n</span><span class="p">)</span> <span class="o">-</span> <span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="p">(</span><span class="n">xmin</span><span class="o">/</span><span class="n">b</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="n">alpha</span><span class="o">-</span><span class="mi">1</span><span class="p">)))</span> <span class="k">for</span> <span class="n">i</span><span class="p">,</span><span class="n">b</span> <span class="ow">in</span> <span class="nb">zip</span><span class="p">(</span><span class="nb">xrange</span><span class="p">(</span><span class="n">n</span><span class="p">),</span><span class="n">z</span><span class="p">)]</span> <span class="p">)</span>
        <span class="c"># Parallels Eqn 3.5 in Clauset et al 2009, but zeta(alpha, xmin) = (alpha-1)/xmin.  Really is Eqn B3 in paper.</span>
        <span class="c">#L = n*log((alpha-1)/xmin) - alpha*sum(log(z/xmin))</span>
        <span class="n">sl</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">([</span><span class="n">math</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">a</span><span class="o">/</span><span class="n">xmin</span><span class="p">)</span> <span class="k">for</span> <span class="n">a</span> <span class="ow">in</span> <span class="n">z</span><span class="p">])</span>
        <span class="n">L</span> <span class="o">=</span> <span class="p">(</span><span class="n">n</span><span class="o">*</span><span class="n">math</span><span class="o">.</span><span class="n">log</span><span class="p">((</span><span class="n">alpha</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span><span class="o">/</span><span class="n">xmin</span><span class="p">)</span> <span class="o">-</span> <span class="n">alpha</span><span class="o">*</span><span class="n">sl</span><span class="p">)</span>
        <span class="c">#requires another map... Larr = arange(len(unique(x))) * log((av-1)/unique(x)) - av*sum</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">_likelihood</span> <span class="o">=</span> <span class="n">L</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">_xmin</span> <span class="o">=</span> <span class="n">xmin</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">_xmins</span> <span class="o">=</span> <span class="n">xmins</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">_alpha</span><span class="o">=</span> <span class="n">alpha</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">_alphaerr</span> <span class="o">=</span> <span class="p">(</span><span class="n">alpha</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span><span class="o">/</span><span class="n">math</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">_ks</span> <span class="o">=</span> <span class="n">ks</span>  <span class="c"># this ks statistic may not have the same value as min(dat) because of unique()</span>
        <span class="c">#if scipyOK: self._ks_prob = scipy.stats.kstwobign.sf(ks*numpy.sqrt(n))</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">_ngtx</span> <span class="o">=</span> <span class="n">n</span>
        <span class="k">if</span> <span class="n">math</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">L</span><span class="p">)</span> <span class="ow">or</span> <span class="n">math</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">xmin</span><span class="p">)</span> <span class="ow">or</span> <span class="n">math</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">alpha</span><span class="p">):</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&quot;plfit failed; returned a nan&quot;</span><span class="p">)</span>

        <span class="k">if</span> <span class="ow">not</span> <span class="n">quiet</span><span class="p">:</span>
            <span class="k">if</span> <span class="n">verbose</span><span class="p">:</span> <span class="k">print</span> <span class="s">&quot;The lowest value included in the power-law fit, &quot;</span><span class="p">,</span>
            <span class="k">print</span> <span class="s">&quot;xmin: </span><span class="si">%g</span><span class="s">&quot;</span> <span class="o">%</span> <span class="n">xmin</span><span class="p">,</span>
            <span class="k">if</span> <span class="n">verbose</span><span class="p">:</span> <span class="k">print</span> <span class="s">&quot;</span><span class="se">\n</span><span class="s">The number of values above xmin, &quot;</span><span class="p">,</span>
            <span class="k">print</span> <span class="s">&quot;n(&gt;xmin): </span><span class="si">%i</span><span class="s">&quot;</span> <span class="o">%</span> <span class="n">n</span><span class="p">,</span>
            <span class="k">if</span> <span class="n">verbose</span><span class="p">:</span> <span class="k">print</span> <span class="s">&quot;</span><span class="se">\n</span><span class="s">The derived power-law alpha (p(x)~x^-alpha) with MLE-derived error, &quot;</span><span class="p">,</span>
            <span class="k">print</span> <span class="s">&quot;alpha: </span><span class="si">%g</span><span class="s"> +/- </span><span class="si">%g</span><span class="s">  &quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">alpha</span><span class="p">,</span><span class="bp">self</span><span class="o">.</span><span class="n">_alphaerr</span><span class="p">),</span> 
            <span class="k">if</span> <span class="n">verbose</span><span class="p">:</span> <span class="k">print</span> <span class="s">&quot;</span><span class="se">\n</span><span class="s">The log of the Likelihood (the maximized parameter), &quot;</span><span class="p">,</span>
            <span class="k">print</span> <span class="s">&quot;Log-Likelihood: </span><span class="si">%g</span><span class="s">  &quot;</span> <span class="o">%</span> <span class="n">L</span><span class="p">,</span>
            <span class="k">if</span> <span class="n">verbose</span><span class="p">:</span> <span class="k">print</span> <span class="s">&quot;</span><span class="se">\n</span><span class="s">The KS-test statistic between the best-fit power-law and the data, &quot;</span><span class="p">,</span>
            <span class="k">print</span> <span class="s">&quot;ks: </span><span class="si">%g</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">ks</span><span class="p">)</span>

        <span class="k">return</span> <span class="n">xmin</span><span class="p">,</span><span class="n">alpha</span>

</div></div>
<div class="viewcode-block" id="plexp"><a class="viewcode-back" href="../../plfit.html#plfit.plfit_py.plexp">[docs]</a><span class="k">def</span> <span class="nf">plexp</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">xm</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span><span class="n">a</span><span class="o">=</span><span class="mf">2.5</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    CDF(x) for the piecewise distribution exponential x&lt;xmin, powerlaw x&gt;=xmin</span>
<span class="sd">    This is the CDF version of the distributions drawn in fig 3.4a of Clauset et al.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">C</span> <span class="o">=</span> <span class="mi">1</span><span class="o">/</span><span class="p">(</span><span class="o">-</span><span class="n">xm</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">a</span><span class="p">)</span> <span class="o">-</span> <span class="n">xm</span><span class="o">/</span><span class="n">a</span> <span class="o">+</span> <span class="n">math</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="p">)</span><span class="o">*</span><span class="n">xm</span><span class="o">/</span><span class="n">a</span><span class="p">)</span>
    <span class="n">Ppl</span> <span class="o">=</span> <span class="k">lambda</span><span class="p">(</span><span class="n">X</span><span class="p">):</span> <span class="mi">1</span><span class="o">+</span><span class="n">C</span><span class="o">*</span><span class="p">(</span><span class="n">xm</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">a</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="n">X</span><span class="o">/</span><span class="n">xm</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">a</span><span class="p">))</span>
    <span class="n">Pexp</span> <span class="o">=</span> <span class="k">lambda</span><span class="p">(</span><span class="n">X</span><span class="p">):</span> <span class="n">C</span><span class="o">*</span><span class="n">xm</span><span class="o">/</span><span class="n">a</span><span class="o">*</span><span class="n">math</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="p">)</span><span class="o">-</span><span class="n">C</span><span class="o">*</span><span class="p">(</span><span class="n">xm</span><span class="o">/</span><span class="n">a</span><span class="p">)</span><span class="o">*</span><span class="n">math</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">a</span><span class="o">*</span><span class="p">(</span><span class="n">X</span><span class="o">/</span><span class="n">xm</span><span class="o">-</span><span class="mi">1</span><span class="p">))</span>
    <span class="n">d</span><span class="o">=</span><span class="n">Ppl</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
    <span class="n">d</span><span class="p">[</span><span class="n">x</span><span class="o">&lt;</span><span class="n">xm</span><span class="p">]</span><span class="o">=</span><span class="n">Pexp</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">d</span>
</div>
<div class="viewcode-block" id="plexp_inv"><a class="viewcode-back" href="../../plfit.html#plfit.plfit_py.plexp_inv">[docs]</a><span class="k">def</span> <span class="nf">plexp_inv</span><span class="p">(</span><span class="n">P</span><span class="p">,</span><span class="n">xm</span><span class="p">,</span><span class="n">a</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Inverse CDF for a piecewise PDF as defined in eqn. 3.10</span>
<span class="sd">    of Clauset et al.  </span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">C</span> <span class="o">=</span> <span class="mi">1</span><span class="o">/</span><span class="p">(</span><span class="o">-</span><span class="n">xm</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">a</span><span class="p">)</span> <span class="o">-</span> <span class="n">xm</span><span class="o">/</span><span class="n">a</span> <span class="o">+</span> <span class="n">math</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="p">)</span><span class="o">*</span><span class="n">xm</span><span class="o">/</span><span class="n">a</span><span class="p">)</span>
    <span class="n">Pxm</span> <span class="o">=</span> <span class="mi">1</span><span class="o">+</span><span class="n">C</span><span class="o">*</span><span class="p">(</span><span class="n">xm</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">a</span><span class="p">))</span>
    <span class="n">pp</span> <span class="o">=</span> <span class="n">P</span>
    <span class="n">x</span> <span class="o">=</span> <span class="n">xm</span><span class="o">*</span><span class="p">(</span><span class="n">pp</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">a</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">C</span><span class="o">*</span><span class="n">xm</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="mi">1</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">a</span><span class="p">))</span> <span class="k">if</span> <span class="n">pp</span> <span class="o">&gt;=</span> <span class="n">Pxm</span> <span class="k">else</span> <span class="p">(</span><span class="n">math</span><span class="o">.</span><span class="n">log</span><span class="p">(</span> <span class="p">((</span><span class="n">C</span><span class="o">*</span><span class="n">xm</span><span class="o">/</span><span class="n">a</span><span class="p">)</span><span class="o">*</span><span class="n">math</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="p">)</span><span class="o">-</span><span class="n">pp</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">C</span><span class="o">*</span><span class="n">xm</span><span class="o">/</span><span class="n">a</span><span class="p">))</span> <span class="o">-</span> <span class="n">a</span><span class="p">)</span> <span class="o">*</span> <span class="p">(</span><span class="o">-</span><span class="n">xm</span><span class="o">/</span><span class="n">a</span><span class="p">)</span>
    <span class="c">#x[P&gt;=Pxm] = xm*( (P[P&gt;=Pxm]-1) * (1-a)/(C*xm) )**(1/(1-a)) # powerlaw</span>
    <span class="c">#x[P&lt;Pxm] = (math.log( (C*xm/a*math.exp(a)-P[P&lt;Pxm])/(C*xm/a) ) - a) * (-xm/a) # exp</span>

    <span class="k">return</span> <span class="n">x</span>
</div>
<div class="viewcode-block" id="pl_inv"><a class="viewcode-back" href="../../plfit.html#plfit.plfit_py.pl_inv">[docs]</a><span class="k">def</span> <span class="nf">pl_inv</span><span class="p">(</span><span class="n">P</span><span class="p">,</span><span class="n">xm</span><span class="p">,</span><span class="n">a</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot; </span>
<span class="sd">    Inverse CDF for a pure power-law</span>
<span class="sd">    &quot;&quot;&quot;</span>
    
    <span class="n">x</span> <span class="o">=</span> <span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">P</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="mi">1</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">a</span><span class="p">))</span> <span class="o">*</span> <span class="n">xm</span>
    <span class="k">return</span> <span class="n">x</span>
</div>
<div class="viewcode-block" id="test_fitter"><a class="viewcode-back" href="../../plfit.html#plfit.plfit_py.test_fitter">[docs]</a><span class="k">def</span> <span class="nf">test_fitter</span><span class="p">(</span><span class="n">xmin</span><span class="o">=</span><span class="mf">1.0</span><span class="p">,</span> <span class="n">alpha</span><span class="o">=</span><span class="mf">2.5</span><span class="p">,</span> <span class="n">niter</span><span class="o">=</span><span class="mi">500</span><span class="p">,</span> <span class="n">npts</span><span class="o">=</span><span class="mi">1000</span><span class="p">,</span> <span class="n">invcdf</span><span class="o">=</span><span class="n">plexp_inv</span><span class="p">,</span>
        <span class="n">quiet</span><span class="o">=</span><span class="bp">True</span><span class="p">,</span> <span class="n">silent</span><span class="o">=</span><span class="bp">True</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Tests the power-law fitter </span>

<span class="sd">    Example (fig 3.4b in Clauset et al.):</span>
<span class="sd">    xminin=[0.25,0.5,0.75,1,1.5,2,5,10,50,100]</span>
<span class="sd">    xmarr,af,ksv,nxarr = plfit.test_fitter(xmin=xminin,niter=1,npts=50000)</span>
<span class="sd">    loglog(xminin,xmarr.squeeze(),&#39;x&#39;)</span>

<span class="sd">    Example 2:</span>
<span class="sd">    xminin=[0.25,0.5,0.75,1,1.5,2,5,10,50,100]</span>
<span class="sd">    xmarr,af,ksv,nxarr = plfit.test_fitter(xmin=xminin,niter=10,npts=1000)</span>
<span class="sd">    loglog(xminin,xmarr.mean(axis=0),&#39;x&#39;)</span>

<span class="sd">    Example 3:</span>
<span class="sd">    xmarr,af,ksv,nxarr = plfit.test_fitter(xmin=1.0,niter=1000,npts=1000)</span>
<span class="sd">    hist(xmarr.squeeze());</span>
<span class="sd">    # Test results:</span>
<span class="sd">    # mean(xmarr) = 0.70, median(xmarr)=0.65 std(xmarr)=0.20</span>
<span class="sd">    # mean(af) = 2.51 median(af) = 2.49  std(af)=0.14</span>
<span class="sd">    # biased distribution; far from correct value of xmin but close to correct alpha</span>
<span class="sd">    </span>
<span class="sd">    Example 4:</span>
<span class="sd">    xmarr,af,ksv,nxarr = plfit.test_fitter(xmin=1.0,niter=1000,npts=1000,invcdf=pl_inv)</span>
<span class="sd">    print(&quot;mean(xmarr): %0.2f median(xmarr): %0.2f std(xmarr): %0.2f&quot; % (mean(xmarr),median(xmarr),std(xmarr)))</span>
<span class="sd">    print(&quot;mean(af): %0.2f median(af): %0.2f std(af): %0.2f&quot; % (mean(af),median(af),std(af)))</span>
<span class="sd">    # mean(xmarr): 1.19 median(xmarr): 1.03 std(xmarr): 0.35</span>
<span class="sd">    # mean(af): 2.51 median(af): 2.50 std(af): 0.07</span>

<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">sz</span> <span class="o">=</span> <span class="n">niter</span>
    <span class="n">xmarr</span><span class="p">,</span><span class="n">alphaf_v</span><span class="p">,</span><span class="n">ksv</span><span class="p">,</span><span class="n">nxarr</span> <span class="o">=</span> <span class="p">([</span><span class="mi">0</span><span class="p">]</span><span class="o">*</span><span class="n">sz</span><span class="p">,)</span><span class="o">*</span><span class="mi">4</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">niter</span><span class="p">):</span>
        <span class="n">randarr</span> <span class="o">=</span> <span class="p">[</span><span class="n">random</span><span class="o">.</span><span class="n">random</span><span class="p">()</span> <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">npts</span><span class="p">)]</span>
        <span class="n">fakedata</span> <span class="o">=</span> <span class="p">[</span><span class="n">invcdf</span><span class="p">(</span><span class="n">r</span><span class="p">,</span><span class="n">xmin</span><span class="p">,</span><span class="n">alpha</span><span class="p">)</span> <span class="k">for</span> <span class="n">r</span> <span class="ow">in</span> <span class="n">randarr</span><span class="p">]</span>
        <span class="n">TEST</span> <span class="o">=</span> <span class="n">plfit</span><span class="p">(</span><span class="n">fakedata</span><span class="p">,</span><span class="n">quiet</span><span class="o">=</span><span class="n">quiet</span><span class="p">,</span><span class="n">silent</span><span class="o">=</span><span class="n">silent</span><span class="p">,</span><span class="n">nosmall</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>
        <span class="n">alphaf_v</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">TEST</span><span class="o">.</span><span class="n">_alpha</span>
        <span class="n">ksv</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">TEST</span><span class="o">.</span><span class="n">_ks</span>
        <span class="n">nxarr</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">TEST</span><span class="o">.</span><span class="n">_ngtx</span>
        <span class="n">xmarr</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">TEST</span><span class="o">.</span><span class="n">_xmin</span>

    <span class="k">return</span> <span class="n">xmarr</span><span class="p">,</span><span class="n">alphaf_v</span><span class="p">,</span><span class="n">ksv</span><span class="p">,</span><span class="n">nxarr</span>
</pre></div></div>

          </div>
        </div>
      </div>
        </div>
        <div class="sidebar">
          <h3>Table Of Contents</h3>
          <ul>
<li class="toctree-l1"><a class="reference internal" href="../../agpy.html">Adam Ginsburg&#8217;s Python Code (agpy)</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../image_tools.html">Image Tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../fft_tools.html">AG_fft_tools Package</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../plfit.html">plfit Package</a></li>
</ul>

          <h3 style="margin-top: 1.5em;">Search</h3>
          <form class="search" action="../../search.html" method="get">
            <input type="text" name="q" />
            <input type="submit" value="Go" />
            <input type="hidden" name="check_keywords" value="yes" />
            <input type="hidden" name="area" value="default" />
          </form>
          <p class="searchtip" style="font-size: 90%">
            Enter search terms or a module, class or function name.
          </p>
        </div>
        <div class="clearer"></div>
      </div>
    </div>

    <div class="footer">
      &copy; Copyright 2011, Adam Ginsburg.
      Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 1.2pre.
    <script type="text/javascript">

      var _gaq = _gaq || [];
      _gaq.push(['_setAccount', 'UA-6253248-2']);
      _gaq.push(['_trackPageview']);

      (function() {
        var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true;
        ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js';
        var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s);
      })();

    </script>
        
    </div>
  </body>
</html>